
#include "errlog.h"
#include "grid3d.h"
#include "g3dvec.h"
#include "glovar.h"
#include "vecs.h"
#include "fe.h"
#include <stdio.h>
#include <string.h>

extern int NO_NORMS;
extern const char nrms_fname[];

/* The es_norms_xxx functions are used when the exact solution is known, i.e.
 * in a convergence test.
 * The td_norms_xxx functions are used when the exact solution is not known and
 * we are solving for a stationary solution.
 */

void td_norms_timestep(double *N,
                       unsigned elemU, const vector u, const vector u0,
                       unsigned elemP, const scalar p, const scalar p0)
{
    N[0]=norm_vec_diff_l2(grid, elemU, u, u0)/TStep;
    N[2]=norm_vec_diff_h1(grid, elemU, u, u0)/TStep;
    N[1]=norm_scal_diff_l2(grid, elemP, p, p0)/TStep;
    N[3]=norm_scal_diff_h1(grid, elemP, p, p0)/TStep;
    PrintLogNorm("u-u0", "L2", N[0]);
    PrintLogNorm("p-p0", "L2", N[1]);
    PrintLogNorm("u-u0", "H1", N[2]);
    PrintLogNorm("p-p0", "H1", N[3]);
    PRINTF("||u-u0||_L2 = %10.4e", N[0]);
}

static void es_norms_print(const double *N);

void es_norms_init(double *N)
{
    memset(N, 0, 5*sizeof(double));
}

void es_norms_timestep(double *N,
                       unsigned elemU, const vector u, const vector exact_u,
                       unsigned elemP, const scalar p, const scalar exact_p)
{
    double norms[4];

    norms[0]=norm_vec_diff_l2(grid, elemU, u, exact_u);
    norms[2]=norm_vec_diff_h1(grid, elemU, u, exact_u);
    norms[1]=norm_scal_diff_l2(grid, elemP, p, exact_p);
    norms[3]=norm_scal_diff_h1(grid, elemP, p, exact_p);
    PrintLogNorm("u-eu", "L2", norms[0]);
    PrintLogNorm("p-ep", "L2", norms[1]);
    PrintLogNorm("u-eu", "H1", norms[2]);
    PrintLogNorm("p-ep", "H1", norms[3]);
    PRINTF("||u-eu||_L2 = %10.4e", norms[0]);
    N[0]=fmax(N[0], norms[0]);
    N[1]=fmax(N[1], norms[1]);
    N[2]=fmax(N[2], norms[2]);
    N[3]=fmax(N[3], norms[3]);
}

void es_norms_done(double *N)
{
    /*
        N[0] = sqrt(TStep*N[0]);
        N[1] = sqrt(TStep*N[1]);
        N[2] = sqrt(TStep*N[2]);
        N[3] = sqrt(TStep*N[3]);
     */
    PrintLogLine(1, ' ');
    PrintLogLine(78, '=');
    PrintLogNorm("u-eu", "L2linf", N[0]);
    PrintLogNorm("p-ep", "L2linf", N[1]);
    PrintLogNorm("u-eu", "H1linf", N[2]);
    PrintLogNorm("p-ep", "H1linf", N[3]);
    PrintLogLine(78, '=');
    PRINTF("||u-eu||_L2linf = %10.4e\n", N[0]);
    PRINTF("||p-ep||_L2linf = %10.4e\n", N[1]);
    PRINTF("||u-eu||_H1linf = %10.4e\n", N[2]);
    PRINTF("||p-ep||_H1linf = %10.4e\n", N[3]);

    if(flags.convergence)
        es_norms_print(N);
}

static void es_norms_print(const double *N)
{
    FILE *F=fopen(nrms_fname, "a+t");
    if(F==NULL)
        return;
    fseek(F, 0, SEEK_SET);
    if(fgetc(F)==EOF)
    {
        fprintf(F, "%-12s    %-12s    %-14s  %-14s  %-14s  %-14s\n",
                "dT", "grid", "||u-eu||L2linf", "||p-ep||L2linf",
                "||u-eu||H1linf", "||p-ep||H1linf");
    }
    else
        fseek(F, 0, SEEK_END);
    /* fprintf(F, "%-12g    %-12g", (double) TStep, (double) GStep); */
    fprintf(F, "%-12g    %-12s", (double) TStep, grid->gridname);
    fprintf(F, "    %12e", N[0]);
    fprintf(F, "    %12e", N[1]);
    fprintf(F, "    %12e", N[2]);
    fprintf(F, "    %12e", N[3]);
    fprintf(F, "\n");
    fclose(F);
}



